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i-^ ' Abstract 

O-i 

A wide class of binary-state dynamics on networks — including, for example, the voter 
model, the Bass diffusion model, and threshold models — can be described in terms of transition 
rates (spin-flip probabilities) that depend on the number of nearest neighbors in each of the 
two possible states. High-accuracy approximations for the emergent dynamics of such models 
on uncorrelated, infinite networks are given by recently-developed compartmental models or 
c/3 ■ approximate master equations (AME). Pair approximations (PA) and mean- field theories can 

be systematically derived from the AME. We show that PA and AME solutions can coincide 
under certain circumstances, and numerical simulations confirm that PA is highly accurate in 
these cases. For monotone dynamics (where transitions out of one nodal state are impossible, 
e.g., SI disease-spread or Bass diffusion), PA and AME give identical results for the fraction 
of nodes in the infected (active) state for all time, provided the rate of infection depends 
linearly on the number of infected neighbors. In the more general non-monotone case, we 
derive a condition — that proves equivalent to a detailed balance condition on the dynamics — 
for PA and AME solutions to coincide in the limit t — > oo. This permits bifurcation analysis, 
yielding explicit expressions for the critical (ferromagnetic/paramagnetic transition) point of 
such dynamics, closely analogous to the critical temperature of the Ising spin model. Finally, 
the AME for threshold models of propagation is shown to reduce to just two differential 
equations, and to give excellent agreement with numerical simulations. As part of this work, 
Matlab code for implementing and solving the differential equation systems is made available 
for download. 

1 Introduction 

Binary-state dynamics on complex networks are frequently used as simple models of social inter- 
action [U [H [31 H]. Each node in a network is considered to be in one of two possible states (e.g., 
susceptible or infected, inactive or active) at each moment of time. Switching of nodes to the oppo- 
site state occurs stochastically, with probabilities that depend on the state of the updating node, 
and the states of its neighbors. Examples include models for competition between two opinions in 
a network-based population, such as the voter model 0[6] and the majority- vote model [SJ [71 E] • 
Models for rumor propagation, adoption of new technologies, or the spreading of behavior are 
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often based on generalizations of disease-spread [lOj [TT] or individual-threshold |T2l [131 03] 
dynamics. Recently, data from online social networks has been mined to help increase the un- 
derstanding of how individuals are influenced by their network neighbors [HI [T71 [181 021 USB I2T] : 
further insight is provided by experiments where the network topology and social interactions 
are carefully controlled [22] . Binary-state dynamics provide the simplest test for newly-developed 
models [23, 24J; interest is often focussed on the critical parameters where the behavior of the 
model changes qualitatively, as at the paramagnetic/ferromagnetic phase transition that occurs 
in the Ising spin model at the critical temperature [251 [2S1 EZ] • 

All these models, and many others, can be considered as special cases of a general formulation 
for stochastic binary-state dynamics. Letting k be the degree (number of network neighbors) of a 
node X in an unweighted, undirected, uncorrelated, static network, we denote by m the number of 
neighbors of X who are in state 1, i.e., infected/active/spin-up, depending on the model context. 
If node X itself is in state (i.e., susceptible/inactive/spin-down) then the rate F^^ m is defined by 
letting Fk iin dt be the probability that X will switch from state to state 1 in an infinitesimally 
small time interval dt. On the other hand, if node X is in state 1, then R^^dt is the probability 
that node X will switch to state within the time interval dt. We call the functions Fk m and Rk, m 
the infection rate and recovery rate respectively; collectively they are also called the transition 
rates or spin-flip rates of the dynamics. These rates depend (only) on the degree k of node X, 
and the number m of neighbors of X who are in state 1; we also assume that the rates m and 
Rk,m do not vary in time. In Sec. [2] (see Table H]) we show that many stochastic binary-state 
social interaction models can be expressed in terms of suitable Fk m and Rk, m functions [281 129j. 
Deterministic threshold models |12 [I131 [T4] — where nodes change state (become active) only when 
the number m of active neighbors exceeds a pre-defined threshold level — may also be put into 
this form (see Sec. [6]). 

Approximations for macroscopic quantities of interest — such as the expected fraction of active 
nodes in the network at a given time — are needed in order to understand how the combination 
of network topology (e.g., the degree distribution P^) and the microscopic dynamics (defined by 
the rates Fk^ m and Rk.m) affect the emergent dynamics at macroscopic level. In the limit of 
infinite network size, systems of deterministic differential equations may be derived to describe 
such macroscopic quantities, at various levels of approximation. The most common analytical 
approach for dynamics on complex networks is mean- field (MF) theory. Mean- field theories are 
typically derived under a number of assumptions, the most important of which for the current 
discussion is the assumed lack of dynamical correlation^. Under this assumption, the neighbors of 
a node A are considered to be active (i.e., in state 1) with a probability that is determined by the 
overall fraction of active nodes (of same degree) across the network as a whole. In particular, the 
probability of a neighbor of A being active is assumed to be independent of the state of node A: 
thus there is assumed to be no correlation between the state of A and the state of ^4's neighbors. 

Examples of MF theories for binary-state dynamics are found in [91 l32"l [3] , and the validity 
of the assumptions of mean-field theories on networks are considered in detail in |33(. [34] . Mean- 
field theories, though relatively simple to derive and solve, are often quite inaccurate, especially 
on sparse networks (with low mean degree z) and near to critical points of the dynamics. Pair 

1 Dynamical correlations should not be confused with degree-degree (structural) correlations, which are a property 
of the network rather than of the dynamics. In this paper we assume zero structural correlations in the networks 
(uncorrelated configuration model [301 1311 14]). 
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approximation (PA) theories improve upon the accuracy of MF theories — at the cost of extra 
complexity in the resulting system of differential equations — by including dynamical correlations 
at a pairwise level [351 E3 EH [381 133 US HI]- In addition to the fraction of active nodes in 
each degree class (as in MF), PA theories account for the fraction of edges in the network that 
link active nodes to active nodes (called 1-1 edges here), inactive nodes to inactive nodes (S-S 
edges), and inactive to active nodes (S-I edges). However, dynamical correlations beyond nearest 
neighbors are not captured by PA theories. If node A is inactive, for example, the probability 
that each of its neighbors is in the active state is approximated in PA theory by the relative 
preponderance of S-I edges over S-S edges in the whole network. Moreover, each neighbor of A is 
considered to be independent of every other neighbor of A. This means that if m is the number 
of A's neighbors that are in the active state, m is assumed to have a binomial distribution. 

Higher-order accuracy, beyond PA level, has been obtained recently using compartmental mod- 
els or approximate master equations (AME) |42l 143} 144" ! 145] . These approaches consider the state 
of nodes and their immediate neighbors, generating large systems of differential equations, which 
better capture the (non-binomial) distributions of the number of active neighbors of updating 
nodes. The increased complexity of the differential equations systems gives improved accuracy 
over PA and MF, particularly near critical points of the dynamics [44] . 

In this paper, we investigate the possibility of, for certain classes of dynamics, reducing the 
dimensionality of the full AME system to obtain a simpler set of differential equations, but without 
any loss of accuracy. Ideally, the reduced-dimension system would permit analytical solutions, but 
even if closed-form solutions are not possible, the smaller system is less computationally expensive 
to solve. We restrict our attention to static networks, but note that AME approaches have also 
been successfully used for models where the network structure co-evolves with the dynamics 
[42 \ l46l 147] . We also assume here that the networks are generated by the configuration model 
[30 } 1311 H]. so having zero clustering (transitivity) in the infinite-size limit. These simplifications 
enable us to focus on the types of dynamics for which the AME can be exactly reduced to a 
simpler (e.g., pair approximation) system; we anticipate that complicating factors such as non- 
zero clustering can be studied in subsequent work, for example by employing the models of 

[MlllSlEolEI]. 

The remainder of the paper is organized as follows. In Sec. [2] we give examples of binary- 
state dynamics and in Sec. [3] we briefly review the approximate master equation approach of [44J, 
showing how the equations of PA and MF theories can be derived as systematic approximations 
of the full AME. In Sec. [4] the AME and PA solutions are shown to give identical results for 
certain macroscopic quantities within the class of dynamics we call monotone. In Sec. [5] a more 
general class of dynamics (corresponding to equilibrium spin systems obeying detailed balance) 
is shown to have AME and PA solutions that are equal in the steady-state limit t — > oo, but not 
at finite t. Focussing on spin systems with up-down symmetry, this result enables us to derive 
an analytical result for the critical point of such dynamics, i.e., the bifurcation point marking the 
phase transition from paramagnetic (disordered) to ferromagnetic (ordered) behavior, as in the 
Ising model on complex networks |26|. 125] . Finally, in Sec. El we show that (monotone) threshold 
models may be accurately approximated using an extended version of the AME, and that the 
large system of equations may be reduced to just two differential equations to determine the 
expected fraction of active nodes. Details of some mathematical derivations are contained in the 
Appendices. 
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Figure 1: Schematic of the infection rate Fk t m an d recovery rate Rkm as defined in the text. 
Here — and in Figs. 01 and flUI — white nodes are susceptible/inactive/spin-down and black nodes 
are infected/active/spin-up. Examples of the dependence of the rates on the number m of infected 
neighbors: (a) susceptible-infected-susceptible disease spread model, (b) Glauber dynamics for the 
Ising model, (c) a stylized "friend-saturation" model. Case (c) is motivated by the data analysis of 
social contagion for Digg stories in Fig. 4 of [T7], where the probability of, for example, becoming 
infected (voting on a Digg story) initially increases with the number of infected neighbors, but then 
saturates, and eventually decreases, giving a function that is roughly symmetric about m = k/2. 



2 Examples of binary-state dynamics 

Here we briefly examine some of the binary-state models that can be described using transition 
rates i*fc m and Rk,m] see Table Q] and the schematic Fig. Q3 Note that in many of the models 
described here a single randomly-chosen node is updated at each time step, for which the time 
increment is dt = 1 /N, where N is the number of nodes in the network (and we take the N — > oo 
limit). In such a model, Fp. m is the probability that a spin-down node, when selected for updating, 
flips to spin-up, while Rk, m is the corresponding probability of a node flipping from spin-up to 
spin-dowrJl. 

The susceptible-infected-susceptible (SIS) disease spread model [52j [531 G3] assumes that in- 
fected individuals transmit disease to each of their network neighbors at a rate A. Thus if a 
susceptible node has m infected neighbors, the probability of it remaining susceptible for a time 



2 The case of symmetric spin models, where the dynamics are unchanged by simultaneous flipping of all nodes' 
spins, is considered in detail in Sec. 15.11 
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Table 1: Infection and recovery rates for some examples of binary-state dynamics on networks: k 
is the node's degree, m is its number of infected neighbors, z is the mean degree (k) = ^2 k kPt 
of the network. The parameters of the various models are described in Sec. [2j Symmetric spin 
models are those with rates that obey condition (f23|) ; equilibrium models obey condition (fT7|) . 
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interval dt is (1 — Xdt) m . Therefore its probability of infection is 1 — (1 — Xdt) m , and in the limit 
dt — > this is Xmdt, giving the Fk^m rate for SIS dynamics in Table [TJ The recovery rate in the 
SIS model is a constant /i, and in the case \x = recovery is impossible; this special case is called 
the SI model. 

The Bass model [55J for diffusion of technological innovations is an extension of the SI model, 
and may be similarly specialized to a network context. In the Bass model, nodes move from the 
susceptible (non-adopted) state to the infected (adopted) state with rate Fk, m = c + dm. The 
parameter c represents individual action, independent of the influence of neighbors, while the 
parameters d quantifies the "herding behavior", whereby individuals copy the actions of their 
neighbors: in this case, by adopting the innovation when m of their neighbors have already done 
so. Transitions from the adopted state to the non-adopted state are not permitted, so Rk,m = 0. 
Kirman's ant colony model |56|, [5TJ EE] for choice between stock market trading strategies has 
similar herding effects, but with switching permitted between both states (strategies). A node 
with m of its k neighbors infected has k — m neighbors in the susceptible state, so the rate for 
transitions of such a node from infected to susceptible is Rk, m = C2 + d(k — m). Note that 
mean-field theory (see Eq. ([5]) below) for the complete-graph case = S^^n — where all nodes are 
connected to all other nodes — gives the standard population-level version of the Bass and Kirman 
models [55], [56] . 

In the voter model [5j [6], a node is chosen at random in each time step (dt = 1/N) and it 
adopts the same state of one of its (randomly-chosen) neighbors. If m of the node's k neighbors 
are in the active state, the node thus becomes active with probability m/k, and becomes inactive 
with probability (k — m)/k, these being the fractions of its neighbors in the respective states. 
Many variants of the voter model have been studied [U [3]. In the link update voter model, for 
example, an edge (rather than a node) is chosen at random in each time step, and one of the 
end-nodes of the chosen edge then copies the state of the other. Nonlinear voter models and 
other variants have been examined by a number of authors |59[ IBTJl l29l I41j . In [41] for example, 
the current state of the node is incorporated into the probabilities for switching state, lending 
an effective inertia to the dynamics: see the transition rates for this case (on 4-regular random 
graphs) in Table [TJ The language model of [60] has been examined in detail in [TO]; in this model 
the two states represent the primary language choice of a person (node), and the probability of 
switching states is proportional to the fraction of speakers in the locality, raised to the power a, 
multiplied by the status parameter s or 1 — s of the respective language (s = 1/2 for the symmetric 
case of equal-status languages, a = 1 then gives the voter model). 

Many models of opinion dynamics are based on the classical Ising model of magnetic spins. 
Here each node is in either the spin-up or spin-down state and transitions occur according to 
dynamical rules which minimise the Hamiltonian of ferromagnetic interactions [3j [61] . Letting 
T represent the temperature of the heat bath, and J the ferromagnetic interaction parameter, 
the transition rates for Glauber [62] and Metropolis [63] dynamics are shown in Table [TJ The 
majority-vote model [5j [7j is an non-equilibrium spin model, with spins tending to align with the 
local neighborhood majority, but with a noise parameter Q giving the probability of misalignment. 

As a final example, we consider threshold models \12\ [T~3"l [HI 164] , which are used to model 
diffusion of fads, collective active, or adoption of innovations [UJ. Each node has a (frozen) 
threshold level, which may depend on the degree of the node. In each time step, a fraction dt 
of the N nodes in the network are updated. When updated, a node compares the number m 
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of its active neighbors to its threshold and activates (with probability 1) if m is greater than or 
equal to its threshold. Similar to the Bass model of innovation diffusion, usually no recovery is 
permitted in this model [T8] (but see [15] for a recent extension including recovery). Note 
that a coordination game, modelling the diffusion of new behavior through a network, can also 
be written as a threshold model of this type |65 } I18 |. [66] . Threshold models, and the extension of 
the AME required to describe them, are considered in Sec. [6] below. 



3 Approximate master equations, pair approximations, and mean- 
field theory 

For completeness, we here briefly recapitulate the approximate master equations (AME) intro- 
duced in [44] : these were derived by generalizing the approach used in |4"3] for SIS dynamics, 
see Appendix. lAl for details. Let s k:m {t) (resp. i k , m (t)) be the fraction of fc-degree nodes that 
are susceptible (resp. infected) at time t, and have m infected neighbors. Then the fraction pk(t) 
of fc-degree nodes that are infected at time t is given by p k (t) = Ylm=o^,m = 1 — ^=0^,™! 
and the fraction of infected nodes in the whole network is found by summing over all £;-classes: 
p(t) = {pk(t)) = Ylk Pk(t). (Recall P k is the degree distribution, i.e., the probability that a 
randomly-chosen node has k neighbors). 

The approximate master equations for the evolution of s k>rn (t) and ik,m(t) are [44] : 

-Jj. S k,m = -F k ,mS k ,m + Rk, m ik,m ~ ft* (k ~ m)s k>m + (3 S (k - TO + l)sfc )m _l 

- 7 s ms ktm + 7 s (m + l)s k , m +i, (1) 
^fc,m = -Rk, m ik,m + Fk, m s k , m ~ P l {k - m)i k ^ m + fi l (k -m + l)i k>m -i 

-l % mi Km + f{m+l)i Km+l , (2) 



for each m in the range 0, . . . , k, and for each £;-class in the network, with (3 £ — — j— 



(Em=o( fc ~" 1 )- F fc,m S k ,m) 
(Em=o( fc - m ) S fe,m) 



s = <E^ (fc-m)H., m i fcm ) < = (J2t omF KmSk ) ^ = (E^uM^ Equations © 

and ([2]), with the time-dependent rates [3 s , 7 s , /3 l and 7* (defined as nonlinear functions of s k ^ m 
and i km ), form a closed system of deterministic equations which can be solved numerically us- 
ing standard methods; Matlab scripts are available from the author's webpage [67] • Assum- 
ing a randomly-chosen fraction ^(0) of nodes are initially infected, the initial conditions are 
Sfe, m (0) = (1 - P(0)) B k , m {p(0)), i k ,m{0) = p(0)B kj7n (p(0)), where B k>m (q) denotes the binomial 
factor 

B k , m (q)=Q^q m (l-q) k - m . (3) 

For dynamics on a general network, with non-empty degree classes from k = up to a cutoff 
kmux, the number of differential equations in the system ([I])-© is (k max + l)(k max + 2), and so 
grows with the square of the largest degree. Some approximation is therefore necessary if it is 
desirable to reduce the AME to a lower-dimensional system. One possibility is to consider the 
parameters p k (t) (resp. q k (t)), defined as the probability that a randomly-chosen neighbor of a 
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susceptible (resp. infected) fc-degree node is infected at time t. Noting that Pk(t) can be expressed 
in terms of s^.m as X^m=o m ^k,m /Y^,m=o ks k,m, an evolution equation for p k may be derived by 
multiplying equation ([I]) by m and summing over m. The right-hand-side of the resulting equation 
contains higher moments of Sk t m, so a closure approximation is needed to proceed. If we make the 
ansatz that s kt m an d i k ,m are proportional to binomial distributions: s k ,m ~ (1 — p k ) -Bfc,m(Pfc)) 
ik,m ~ Pk Bk t m(qk),, we obtain the pair approximation (PA), consisting of the 3/c max + l differential 
equations: 



d 

dt pk = 


k 

-pk Y 


Rk,m 


d 

Jt Pk = 


k 

Y h 

m=0 


m~ 
~ "fe- 


d 

dt qk = 


k 

£ [« - 


rn" 



^ Rk,mBk,m(Qk) + 0- ~ Pk) Y/ F k, m Bk,m{Pk) , 



m=0 



Fk,mBk,m{Pk) 
Rk,mBk,m{qk) 



m=0 
Pk 



Pk 
Pk 



Pk 



Rk,mBk : m{qk j 
Fk,m,Bk,m(Pk) 



+ P S (l~Pk) 

+ F0--qk) 



7 Pfc, 



7 ?fe 5 



(4) 



for each £;-class. The rates here are given by inserting the binomial ansatz into the general 
formulas, so that /3 s , for example, is {(1 - p k )^ m (k - m)F k>m B k!m (j) k )) / {(1 - p k )k(l - p k ))\ 
initial conditions are pfc(O) = p k (0) = 0^(0) = p(0). 

A cruder, mean-field (MF), approximation results from replacing both pk and qk with u: 
« (1 - /9 fc ) B k>m (uj), Bk,m(u), where w = (fpfe) is the probability that one end of 

a randomly-chosen edge is infected, and z = (k) is the mean degree of the network. Using this 
ansatz in the master equations yields a closed system of fc max + 1 differential equations for the 
fraction pk of infected fc-degree nodes: 



d 



-^Pk — —Pk Yj B-k,mBk,m{u) + (1 — p k ) ^ -Ffc,m-Bfc,m(w) 



(5) 



m=0 



m=0 



with P k(0) = p(0). 

In |44| we showed that the AME system (|3ip generally gives improved accuracy over the 
approximations (j4|) and ([5]). Moreover, the general equations (|4|) and (JS]) for pair approximation 
and mean-field theory reduce to previously-studied equations when specific rates Fk t m and Rk, m 
are selected from Table [TJ In this paper we concentrate on finding dynamics for which the AME 
system can be reduced exactly to lower-dimensional systems, for instance to the PA equations. 



4 Monotone dynamics 

We first consider the case where Rk, m = for all k and m, with F km nonzero for at least 
some arguments. In this case, it is impossible for nodes to switch from the infected state to the 
susceptible state, and so we call this type of dynamics monoton® In this case the Sk, m equations 

3 The case where Fk >m — with some Rk,m nonzero is obviously equivalent, up to swapping of state labels, to 
the case considered here, and so is not considered separately. 
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of the AME system are decoupled from the i kjTn equations, with Eq. ([T]) reducing to 

^Sfc,m = -Fk,mSk,m ~ '(k - m)s k>m + (3 S (k - 171 + l)sfc, m -l, (6) 

and the fraction of infected nodes is given by p(t) = Yl k P k Ylm=o U ~~ s k,m(t)). Thus several 
important macroscopic quantities (though not all) can be found by solving Eq. © for s km (t), 
without considering the i kjm (t) variables. 

Now consider whether or not the pair approximation ansatz s km (t) = (1 — p k {t))B k ^ m {p k {t)) — 
that was used to derive Eqs. @ — could be an exact solution of Eq. ©. Inserting the PA ansatz 
into Eq. (jBJ) and dividing by (1 — Pk)Bk. m (Pk) gives the condition 

1 dp k ( m k-m\dp k 1 - p k 

+ ~, -77- = -Fk,m ~ P [k-m)+ p m, (7) 



d „ / \ f m k — m\ 
7 Bt, m (p) = - - -, B k>m (p) (8) 



l-pk dt \p k 1-PkJ dt " ' p fc 

where the identities 

dp ^'" iyiJ/ yp i _ p 
and 

(k-m + i)B k ,m(p) = —^-mB k ^ m (p) (9) 

have been respectively utilized on the left-hand side and right-hand side of Eq. ([7]). 

Equation ([7|) can be viewed as a condition on the forms of the infection rate F k m for which 
the PA ansatz on s ktTn is an exact solution of the corresponding approximate master equation. 
Note that all terms in Eq. (|7j) are linear in m, so F km is necessarily of the form 

F k ,m = c fc + d k m, (10) 

for (possibly /c-dependent) constants c k and d k . Using this form in (j4j) confirms that the solutions 
of the AME for s km are given by the PA ansatz, with p k and p k being the solutions of the reduced 
system 

^-Pk = (1 - Pk)(ck + kp k d k ) 
dt 

—Pk = -dkPk{l-Pk) + p S (l-Pk), (11) 

where /3 s = (Y, k Pk{^ ~ Pk)K l ~ Pk)(c k + (k — l)p k d k )) / (J2 k p k0- ~ Pk)k(l ~ Pk))- For the spe- 
cial case of z-regular random graphs (i.e., P k = S k)Z for integer z), system (fTT]l can be solved 
analytically to give the explicit solution for the infected fraction as 



pit) = 1 - (1 - p(0))e' ic+zd)t 



_ (1 - p{0))d(z - 2) / _ {c+d(z _ 2))t \ 
c + d(z-2) \ ) 



(12) 



In Figure [2] we show results for SI dynamics (see Table [T]), which are monotone, and of form 
TU|) with c k = and d k = A for all k. Note the exact match between the AME and PA solutions 
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Figure 2: Susceptible-infected (SI) dynamics (with transmission rate A = 1) on a truncated scale- 
free network, with degree distribution given by Pk oc k~ 2 ' 5 for degrees in the range 3 < k < 20, 
with P k = otherwise. Here, and in all subsequent Figures, symbols show the results of numerical 
simulations on networks of size N = 10 5 , using time step dt = 10 -4 . Results are averages over 24 
realizations; error bars indicate mean ± one standard deviation (but note error bars are smaller 
than the symbols in this example). 

for p(t) and for the fraction of S-I edges (see Eqs. (f39l) and (i56j) ) in the network, and the excellent 
agreement with numerical simulation on networks of size N = 10 5 . 

The parameters dk in (|10p may be negative (provided that all F^ im values are non- negative) . 
As an example, consider the Bass diffusion model on a z-regular random graph with c = 1 and 
d = —1/z. This might model, for example, the adoption by indie music fans of a new band, 
where the attractiveness of the band decreases as the number of neighbors m who have already 
"jumped on the bandwagon" increases; see [H] for another approach to this aspect of social 
contagion modelling, which they call "limited imitation contagion" . The expected fraction of fans 
p(t) is given explicitly by Eq. f)12[) ; note for these parameters the entire network does not become 
infected, see Figure El 

It is interesting to note that while, in the monotone case (|10p . the AME solutions for s/% )m (i) 
are exactly reproduced by the PA equations, the corresponding ik, m (t) variables are not necessarily 
equal to their respective PA ansatz Pk(t)Bk !m (qk(t)) ■ Focussing on the SI model, the effect of the 
non-binomial m distribution in ik,m is n °t visible in Figure O as the quantities shown there can 
be written in terms of Sfc jm only; indeed, it is necessary to examine connected triples of nodes to 
demonstrate this effect [38] . Figure [4] shows the fraction of connected triples (defined by choosing 
a node at random, then randomly choosing two of its neighbors) that are of type S-I-S: the 
chosen (middle) node being infected, while both chosen neighbors are susceptible. In the AME 
formulation, this fraction is given by 2^ fc Pj~ X^m=o(^ ~~ m){k — m — l)ik,m/ {k(k — 1)), whereas 
the corresponding PA version is 2 ^ fc P^p^kik — 1)(1 — qk) 2 / {k(k — 1)). The differences seen 
in Fig. 0] — and the close fit of AME results to numerical simulations — indicate that the match 
between AME and PA results in Figs. [2] and [3] is due to the binomial m distribution Sfc jm for 
susceptible nodes, but the corresponding distribution i^-m for infected nodes is not binomial in 
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t t 

Figure 3: Bass diffusion dynamics on a 4-regular random graph (P*. = 6^4), with parameters 
c = 1 and d = —1/4. The AME and PA solutions for the fraction of adopted nodes are both given 
by the explicit formula ([79]) . 



m. 



5 General dynamics and spin models 

We consider in this section the more general case of non-monotone dynamics, where none of the 
rates m and Rkm are identically zero. Inserting the PA ansatz for m and ik,m 

into the AME 

equations (P)-© we find, after dividing by (1 — Pk)Bk m (Pk)> t ne condition 



Pk B k , m (q k ) _ j!) i j2) n 

1 — Pk Bk,m(Pk 



F k , m + R k ,m , "\ T m Z[ = 4 lj + ^m, (13) 



where dp and c*^ represent combinations of terms that depend on t and k, but are independent 
of m. Similarly, using the PA ansatz in the AME equation for % km yields 

^l^f4rH 3, +4V (14) 

Pk £>k,m{Qk) 

where, like Eq. (I13p . the left-hand side is an exponential function of m (because of the presence of 
the Bk t m functions, see Eq. (J3|)), while the right-hand side is linear in m. Bearing in mind that the 
transition rates Fk, m and Rk, m are time-independent, it is not generally possible to simultaneously 
solve equations f)13|) and f|14[) to obtain constant transition rates. We conclude that in this general 
case, AME solutions are not exactly equal, for all times, to the corresponding PA solutions. 

However, an important special case occurs if we restrict our attention to the steady-state limit 
t — > 00. The AME solutions and PA solutions can be identical in the limit t — > 00 (despite being 
different at finite t) if equations (fT3|) and (fl~4"|) are simultaneously satisfied in the steady state. 
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The m dependence of the left- and right-hand sides then requires that 

(15) 



Fk.m Pk Bk,m{<lk) 



Rk,m I - p k B^miPkY 

and cjj. 1 ^ = = = cjj^ = 0, where the overbar denotes the steady-state limit of the 
corresponding variable, e.g., p k = Hindoo Pk(t)- 

The conditions ci = ci = can be shown to imply that (1 — p k ) = j3 /7 s , and the 
^-independence of the right-hand side implies that p k must, in fact, be independent of k: 

Pk =P for all k. (16) 

Similarly, an analysis of the conditions cj^ = cj^ = shows that the steady-state values of q k 
must also be independent of k, i.e., q k = q for all k. 

Replacing p k and q k in Eq. (|15p with these simplifications, we can rewrite the condition on 
the rates as 

F '"" ■b k a m , (17) 



Rk,: 



where 



and 



- /1 — \ k 

Pk (l-q 



1 - Pk V 1 - v 

q(l -p) 



(18) 



(19) 



p(l - q) 

It is demonstrated in Appendix [B] that Eq. (|17p is precisely the condition for microscopic re- 
versibility (detailed balance) of the stochastic dynamics, i.e., for the dynamics to correspond to 
an equilibrium spin-flip model. Moreover, we show in AppendixOthat the steady solution of the 
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AME for transition rates obeying condition ()1T|) can be fully specified in a simple manner, by the 
equation 

Pk = z _ _ fc , (20) 

h + (l-p + pa) 

where p is a solution of the algebraic equation 

p(l — p + pa) \ - 



l-p 2 +p 2 a z b k + (l-p + pa) 

As an example of dynamics obeying condition (fT7|) . we introduce a modified model of susceptible- 
infected-susceptible type, with 

F k ,m = X m and R k , m = H (22) 

for positive constants A and /i. Like the standard SIS disease-spread model (see Table [1]), recovery 
of infected nodes occurs at a constant rate /i. However, in constrast to the linear-in-m dependence 
of the SIS infection rate i*fc, m , here the probability of infection is assumed to grow exponentially 
with m. While admittedly artificial, such dynamics might find application in fitting data from 
social contagion experiments on networks, where the dependence of the rates of m is not yet clear 
[TBI EU E2]. Figure [5] shows that the AME solutions and PA solutions indeed agree as t — > oo, 
despite being different for finite t. 

Other important examples of equilibrium models obeying (|17() are given by the Glauber and 

4 j 

Metropolis dynamics for the Ising spin model (see Tabled] and Fig. [6j), which both have o = e^r 
and 6fc = e~^r k in Eq. fjlTf) . These models are examples of spin models with up-down symmetry, 
and are considered in greater detail in Sec. 15.11 below. 



13 




Figure 6: Glauber dynamics for the Ising spin model dynamics on a Poisson (Erdos-Renyi) ran- 
dom graph with mean degree z = 7. The interaction parameter J is set to 1, and the temperature 
T is 2/log(2.5) « 2.18. 

5.1 Spin models with up-down symmetry 

Spin models with up-down symmetry have dynamics that are invariant to the swap of state labels 
(susceptible to infected, and vice versa) for all nodes. In [281 [29] this symmetry is called "Z2- 
symmetry"; it is characteristic of the voter model and other opinion dynamics model; see the 
fourth column of Table [TJ In terms of the transition rates, the symmetry condition implies that 

Rk m = Fk,k-m for m = 0, . . . , k and for all k. (23) 

Note that for such dynamics, the AME system ([I])-© is invariant under the change of variables 
Sk,m ^ ik.k-m and ik,m ^ Sfe,fe-m- Since the expected fraction of degree-/c nodes can be written as 
Pk = Ylm=o ^k,m = Sm=o(l ~~ s k,m)i this symmetry condition implies that a solution exists of the 
AME with pk(t) = 1/2 for all t. However, this may not be the only possible solution: depending 
on the initial condition p(0), and on the parameter regime, other solutions of the AME may also 
be found. 

Focussing now on equilibrium spin models with up-down symmetry, which obey condition (|17p 
in addition to (123H . we investigate the stability of the symmetric solution with ~p k = 1/2. First, 
note that there is only a single parameter in these models: putting m = k/2 for even k, or, for odd 
k, m = (k — l)/2 and then m = (k + l)/2, into Eq. (fT7|) and imposing condition (|23|) immediately 
yields the necessary relation 

b k = a-~2 (24) 

between the parameters of equilibrium spin models. Using the steady-state solution of Sec. [5l it 
is possible to show that a critical value a c of the parameter o exists: in the language of dynamical 
systems, this is a (pitchfork) bifurcation point [68J. For parameter values a with a < a c , the 
symmetric solution p = 1/2 is stable, meaning that if p(0) is close to 1/2, the steady-state solution 
will be ~p = 1/2. In spin models (where the magnetization can be written as M = \2p — 1|) this 
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regime is the paramagnetic (disordered) phase; for opinion models in this regime the two opinions 
coexist equally on the network. However, if the parameter a exceeds the critical value a c , then 
the symmetric solution p = 1/2 is unstable, and two other stable solutions, symmetric about 
p = 1/2, exist. This is the ferromagnetic (ordered) phase; for opinion models, one of the two 
opinions dominates the other. The critical value a c gives the phase transition point, and the 
behavior of p near a c can be determined from Eq. (|2ip (see Appendix [D] for details) in a very 
similar fashion to the analysis of the Ising model in [261 125] . see also [69]. The results of such 
analysis (see Appendix [D]) may be summarized as follows. If the degree distribution Pj. possesses 
a finite fourth moment (k ) = k 4 Pk, then the phase transition is of mean-field type, with 
critical parameter 

' ^ V (25) 



c \{k 2 )-2{k)J 

and with p — 1/2 ~ ±(a — a c )2 as a — > a c from above. Following (26j [25j (see Appendix iDj) . 

if the network has a scale- free degree distribution P^ ~ k~^ as k — > oo, then for exponents 7 

1 

in the range 2 < 7 < 3, the critical point is a c = 1, with p — 1/2 ~ ±(a — a c ) 3 ~T, while for 

exponents with 3 < 7 < 5, we have a c given again by Eq. (|25[) . but with near-criticality scaling 

1 

of p — 1/2 ~ ±(a — a c ) 1- 3 as a —7- a c . The case 7 = 3 shows an infinite order transition at a = 1, 
as discussed in [25] . 

As mentioned at the end of Sec. the Ising model (Glauber or Metropolis dynamics) is of 
type ([23]) . with temperature T related to the parameter a via T = 4J/lna, so a = 1 corresponds 
to infinite temperature. A social influence model of this type might be given, for example, by 
transition rates with the properties 

Rk,m = Fk,k-m = Fk,m for m = 0, . . . , fc and all k. (26) 

Here all Fk, m values for m = to |_f J are free parameters; for example, the rates Fk, m an d Rk,m 
could be given by Fig. QJc), which is motivated by the data analysis in Figure 4 of [17J. In any 
model satisfying ([26]) the parameter a of ([23]) is equal to 1, and by the results above the model is 
in the paramagnetic (disordered) phase for all network topologies. However, Eq. (I25p shows that 
a = 1 is near the critical point of the system if the degrees in the network are very heterogeneous 
(so that (fc 2 ) 3> (k)), and indeed the model is poised precisely at criticality (a c = 1) on scale-free 
networks with infinite- variance degree distributions. 

For non-equilibrium spin models, or models not obeying (|17p . the PA and AME solutions are 
different, even as t — > oo, see Fig. [7] for an example using majority-vote dynamics. The solutions 
of the PA equation are therefore of limited usefulness, but it is nevertheless interesting that for 
z-regular graphs an expression for the critical point can be explicitly obtained. In Appendix lEl 
we show that the PA critical point of spin models on such networks occurs precisely when the 
steady-state PA paramater p has the value 

and the location of the critical point in parameter space is given by the solutions F zm (for 
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0.5 



AME 




Figure 7: Majority-vote model on 3-regular random graph, with noise parameter Q = 0.07. This 
non-equilibrium model does not obey condition (I23p and the AME and PA solutions are not equal, 
even as t — > oo, in contrast to Figs. [5]and[6l 



m = 0, . . . , z) of the implicit equation 



m=0 



2m 

z 



F B 



z,m A -'z/m 



2z 



0. 



(28) 



Using the infection rate F^^ m for the majority- vote model (see Table [I]), for example, in Eq. (|28p . 
gives an explicit expression for the PA critical noise parameter Q: 



Qc 



Ylm=0 (l 



2m\ d / z-2 
-) &z,m I 2z-2 



2m 



)B Z 



z-2 
2z-2 



(29) 



Equation (|28p can similarly be used to obtain analytical expressions for the PA critical points in 
other models on z-regular random graphs, e.g. Fig. 12 of [41] . Fig. 1 of |37j . or Fig. 1 of [70] . 



6 Threshold dynamics 

Threshold models are often used to model propagation of fads or collective action through a 
population [T21 [THJ HH HS1 ED IM] • Each node has a (frozen) threshold level; these thresholds may 
be chosen at random (e.g., from a Gaussian distribution), or assigned in some other way (perhaps 
depending on the degree of the node). In the asynchronous-updating version of these models, 
a fraction dt of nodes are randomly chosen for updating in each time step@. When chosen for 
updating, an inactive node becomes active (with probability 1) if m, the number of its neighbors 

4 The extreme case dt = 1 is called synchronous updating |72U15ll73] : however, our interest is in the asynchronous- 
updating case with dt infinitesimally small, which gives continuous-time dynamics in the limit dt — > 0. 
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who are active, exceeds the node's threshold^. Once activated, nodes cannot subsequently become 
inactive, so the dynamics are monotone, cf. Sec. [H Unlike in Eq. (jlOp of Sec. H] however, the 
transition rate is not linear in m; in fact, it is given by 

[1 if m > Mk, 

to reflect the deterministic activation of a node (once it is chosen for updating) when m exceeds 
the threshold level Mk. We have introduced here the vector k to encode two properties defining a 
class of node: their degree k (a scalar), and their type r, which together determine the threshold 
Mk for such nodes. The types are assumed to be from a discrete set of possibilities: all nodes of 
type r = 1, for example, might have the same threshold Mi, with all nodes of type r = 2 having 
a common threshold M2, with M2 7^ M\. In this way, the set of all nodes may be partitioned 
into disjoint sets labelled by their degree and their type; in mathematical notation we combine 
these into a 2-vector, defining k = {k,r} for the k-class of nodes. All nodes in the same k-class 
have the same degree and the same type, and therefore all share the same threshold Mk- We 
generalize the degree distribution P^ to the distribution Pk, which gives the probability that a 
randomly-chosen node has vector k (i.e., has degree k and type r). For example, if the thresholds 
of the nodes are randomly chosen, independent of their degrees, then the Pk distribution can be 
written as P k = PfcPr ; where Pt is the degree distribution and P r is the probability that a node is 
of type r. By taking the discrete set of types to be sufficiently large, it is possible to approximate a 
continuous distribution of types or thresholds with a desired level of accuracy. With this extended 
notation, the AME approach can be generalized in an obvious manner, essentially replacing the 
scalar degree k by the vector k as appropriate in Eq. ([6]), to yield 

^ s k,m = -Fk,m,Sk,m ~ P S (.k - m)s k ,m + P s {k - m + l)s k ,m-i for m = 0, . . . , k, (31) 

with the rate B s given by B s = ^k^ k ^"i=o( fc m )^ m V _ Note the sums here are over all k-classes, 

i.e., over all degrees k and all types r: ^ k := 

In [73] (see also [71]) it is argued that for no-recovery threshold models of the type ([5U1) - ([3T]) 
described above, the fraction p(t) of active nodes at times t can be found by solving just two 
differential equations: 



d P = K4>) ~ P 



dt 
d 

It 



9{ 



J, with 0(O) = /O (O) = ^P k Pk(O), (32) 



where 

K<t>) = ^P k 



p k (0) + (1 - Pk(0)) Yl B k,m(<t>) 



rn>Mu 



(33) 



5 Chapter 19 of 18 shows how a coordination game played on a network — as a model of direct-benefit effects 
from the diffusion of new behavior [65] — can also be expressed as a threshold model of the type considered here. 



17 



AME 

PA 

MF 

• numerical 



5 10 15 20 25 30 35 

t 



0.6 
0.5 

c 0.4 

o 

u 

CS 

"S 0.3 
aa 
— 
o 

« 0.2 
0.1 







AME 

PA 

MF 

• numerical 







10 15 20 25 30 35 
t 



Figure 8: Threshold model on a 5-regular random graph, with all nodes having the same threshold 
M = 2. The AME solution is identical to the solution of the 2-dimensional system (|32l) . 



and 



^ ' 7 



Pk (0) + (1 - p k (0)) Yl B k-i,- 



m>Af k 



(34) 



Here pk(0) is the fraction of nodes with vector k that are activated (infected) at time t = 0; as in 
[73] . we generalize the usual infected fraction p(0) — that is used elsewhere in this paper — to allow 
for possible dependence on the degree or type of the nodes chosen to "seed" the contagion. 

In Appendix |F] we demonstrate that Eqs. ([30]) and (f3Tj) reduce to (f32"]) - (f34"j) through an exact 
solution of (f3"Tj) given by 



Sk,m(t) = (1-Pk(0))£. 



k,m\ 



for m < Mk. 



(35) 



The distribution of Sk, m for m > Mk is, in general, not of the binomial form (|35p . but it can 
nevertheless be given explicitly as detailed in Appendix[Fl This means that the reduced-dimension 
system (j32|) is not precisely of the PA type Q , but it enables efficient and very accurate solution 
of threshold-dynamics models, see the example in Fig. [8l 



7 Conclusions 

We have pointed out the many stochastic binary-state dynamical systems on networks can be 
described using the transition rates F^^ m and Rk, m introduced in Sec. [0 The main results of 
this paper are the identification of dynamics for which the full approximate master equation 
system (f3~TI) can be reduced, without loss of accuracy, to a lower-dimensional system of equations, 
occasionally even yielding closed-form solutions (such as Eqs. (fT2|) and (j28|) ). We showed in Sec.H] 
that the pair approximation system exactly matches the AME results for nodes in one state (e.g., 
the susceptible state) if the dynamics are monotone (i.e., the recovery rate Rk, m is identically 
zero) and the infection rate is linear in the number of infected neighbors (-Pfe m = Ck + d^m). 
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The classical SI model and the Bass diffusion model are of this type — see Figs. [2] and [3] — and 
for z-regular random graphs Eq. (1121) explicitly gives the expected fraction of infected nodes. 
Interestingly, the PA solution is not exactly equal to the AME for nodes in the infected state, see 
Fig.H 

In Sec. [5] we showed that if both sets of transition rates F/. m and Rk, m are non-zero, then 
the AME solutions cannot be given exactly by the PA ansatz for all time t. However, in the 
special case defined by Eq. ([17]), which corresponds to equilibrium stochastic dynamics, i.e., those 
obeying detailed balance (microscopic reversibility), the AME solutions reduce to the PA solution 
in the limit t — > oo, see Figs. [5] and [6] This property allows us to perform bifurcation analysis for 
spin systems with up-down (Z2) symmetry (i.e., for dynamics obeying condition Eq. (|23p ). giving 
the explicit expression (1251) for the critical point where the steady-state limit of p changes from 
the disordered (paramagnetic) state p = 1/2 to an ordered (ferromagnetic) state with p 7^ 1/2. 
The Glauber and Metropolis dynamics for the Ising spin model are important examples that obey 
both conditions (fTTJ) and ([23]) . and our results reproduce the critical Ising temperature that was 
found using very different methods (replica trick, recursion approach) in [261 [25] . The analysis 
of [26] [25] for the critical behavior on scale-free networks can also be applied here, and gives 
results described in Sec. 15.11 and Appendix [Dj We highlight the fact that the critical point a = 1 
for networks with infinite-variance degree distributions is attainable in plausible social influence 
models of this type, as is the case a < 1, although for the Ising model a = 1 corresponds to 
infinite temperature and the regime a < 1 is unphysical. We also give an explicit expression 
(Eq. (|28p) for the critical point (pitchfork bifurcation point) for pair approximation of general 
spin models — including non-equilibrium dynamics that obey (j23[) but not necessarily (|17p — on 
z-regular random graphs, with the caveat that this PA critical point is not identical to the AME 
critical point except for the class of equilibrium models obeying (117D and (123j) . Finally, in Sec. [6] 
we showed that the AME approach can be extended to include threshold models of adoption 
or fad diffusion, and that the AME system can be reduced to a system of only two differential 
equations, see Eqs. ([32]) and Fig. [8] 

The exact agreement of AME and PA solutions — whether for all time as in Sec. H] or in the 
steady-state as in Sec. [5] — implies that higher-order correlations (beyond nearest-neighbor) are 
correctly captured by PA in the cases we have identified. Indeed, the agreement of the theory 
curves with numerical simulation results in all these cases (for all time in Figs. [2] and [3] and 
as t — > 00 in Figs. [5] and [6]) is excellent. We interpret this as indicating that for the classes of 
dynamics where PA is equal to AME, the results of PA are essentially exact. In contrast, in cases 
where PA and AME are not equal, as in Fig. [Jj for example, it is necessary to solve the full AME 
system to obtain high-accuracy approximations. 
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Figure 9: Schematic of transitions to/from the Sfc )m and I^.m sets, as described in equations (|40p 
through (|48p . For each set, the central (Ego) node is shown along with some of its neighbors: 
white nodes are susceptible/inactive/spin-down and black nodes are infected/active/spin-up. See 
also Fig. 1 of 02]. 

A Derivation of Master equations 

We consider binary-state dynamics on static, undirected, connected networks in the limit of infinite 
network size (i.e., N — > oo, where N is the number of nodes in the network). For convenience, we 
call the two possible node states susceptible and infected, as is common in disease-spread models. 
However, this approach also applies to other binary-state dynamics, such as spin systems [37J, 
where each node may be in the +1 (spin-up) or the —1 (spin-down) state. The networks have 
degree distribution and are generated by the configuration model |30t [3Tj 0] . Dynamics are 
stochastic, and are defined by infection and recovery probabilities Fk, m and Rk,m, which depend 
on the degree A; of a node, and on the current number m of infected neighbors of the node, see 
Sec. El 

We now proceed to derive the approximate master equations for dynamics of this type, closely 
following the approach used in [32j [33] for SIS dynamics. Let Sfc m (resp. Ik, m ) be the set of 
nodes which are susceptible (resp. infected), have degree k, and have m infected neighbors. To 
quantify the size of these sets, define m (t) (resp. ik,m(t)) as the fraction of ^-degree nodes that 
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are susceptible (resp. infected) at time t and have m infected neighbors. Then the fraction pk{t) 
of /c-degree nodes that are infected at time t is given by 

k k 

Pk(t) = ^ ik > m = 1 ~~ ^2 Sk ' m i ( 36 ) 

m=0 m=0 

and the fraction of infected nodes in the whole network is found by summing over all /c-classes: 

p{t) = Y J PkPk{t). (37) 

k 

If a randomly-chosen fraction p(0) of nodes are initially infected, then the initial conditions for 
Sk, m and ik, m are easily seen to be 

Sfc,m(0) = (1 - p(0)) Bfc,m(p(0)), ifc,m(0) = ,0(0)^^(^(0)), (38) 

where B k ^ m {q) is the binomial distribution introduced in Eq. ([3]). Note that we can also calculate 
the number of edges of various types using this formalism. For example, the number of edges in 
the network which join a susceptible node to an infected node (we call these S-I edges for short) 
can be expressed in two equivalent ways: 

k k 
N ^2 p k^2 ms k, m or N^Pk ^2(k-m)i ktm . (39) 

k m=0 k m=0 

The first of these expressions, for example, follows from noting that in a sufficiently large network 
that there are NP k nodes of degree k, of which a fraction s kjTn are susceptible and have m 
infected neighbors. Each such node contributes m edges to the total number of S-I edges. Similar 
expressions may also be given for the number of S-S and I-I edges in the network. We note that 
the equivalence of the two expressions in (|39h is preserved by the evolution equations described 
below. 

Next, we examine how the size of the S km set changes in time. We write the general expression 

Sk,m(t + dt) = Sk, m {t) — W(Sk, m —^ Ik,m) Sk, m dt + W(Ik, m —!' Sk, m ) lk,mdt 

— W(Sk, m — > Sk,m+l) s k,m dt + W{Sk,m-l ~^ Sk,m) «fc,m-l dt 
—W(S}.,m — > Sk,m-l) Sk,m dt + W(Sk, m +l ~^ Sk,m) Sfc,m+1 dt (40) 

to reflect all the transitions whose rate is linear in dt (all other state-transitions are negligible in 
the dt — >• limit), see Fig. [9l Here W{Sk^ m — > Ik, m ) dt, for example, is the probability that a node 
in the Sk, m set at time t moves to the Ik, m set by time t + dt. It is clear from the definitions above 
that 

W(S k>m ^Ik tm ) =F k>m and W{I k , m ^S ktm ) = R kjTn . (41) 

A node moves from the S k)7n -i set to the S k ^ m set if it remains susceptible, while one of its 
susceptible neighbors becomes infected. Note this means that an S-S edge changes to an S-I 
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edge. If we suppose that S-S edges change to S-I edges at a (time-dependent) rate /3 s , we can 
write [f| 

W(S ktm ^S ktm+1 ) = (3 s (k-m) and W{S Km ^ S k , m ) = $ s {k - m + 1), (42) 

since nodes in the S km set have k — m susceptible neighbors, while those in the Skm—x set have 
k — m + 1 susceptible neighbors. To calculate /3 s , we count the number of S-S edges in the network 
at time t, and then count the number of edges which switch from being S-S edges to S-I edges in 
the time interval dt; the probability f3 s dt is given by taking the ratio of the latter to the former, 
i.e., 

ps dt = Efc Pk Em=o( fc ~ ™)Fk,m Sk,m dt 
Efc P fcEm=o( fc - m ) s k,m 
A similar approximation is used to define j s , the rate at which S-I edges change to S-S edges 
due to the recovery of an infected node: 

Y _ Efc Pk Em=o(^ ~ m )Rk,m ik,m 
Eifc- P fcEm=o( fc - m Km 

and we then write 

W(S k ,m^S k , m -i) =7 s m and W(S k , m+ i ->S k , m ) = l s {rn + 1). (45) 

Taking the limit dt — > of equation ([40 p gives the master equation for the evolution of s k m (t) 
(see Fig. EJ): 

^Sfc,m = -Fk,mS k , m +R k , m h,m-l3 s (k-m)s kjm +f3 s (k-m+l)s k:m „ 1 -j s m 

(46) 

where m is in the range 0, . . . , k for each fc-class in the network (and adopting the convention 
s k _i = — 0)- Applying identical arguments, mutatis mutandis, to the set I km , we derive 

the corresponding system of equations for i k>m (t): 

4;*k,m = -Rk, m ik, m + Fk,mS k ,m ~/3 l (k- m)i k m + (3\k -m + l)i k , m -i ~ l l mi k)m + f(m + l)i k . m +i, 
dt 

(47) 

for m = 0, . . . ,k and for each /c-class in the network, with time-dependent rates f3 l and 7* defined 
though s kt7n and i k>m as 

0, _ Efc Pk Em=0 m Fk,m S k ,m y _ Efc Pk Em=0 m Pk,m h,m ,^gs 

Efc Pk Em=0 171 S k,m Efc Pk Em=0 171 ^ k ,m 

The approximate master equations (|46p and (|47p . with the time-dependent rates /3 s , j s , f3 l 
and 7* (defined as nonlinear functions of s k)m and ifc >m ), form a closed system of deterministic 



6 The main approximation here is to assume that the edge-state transition rate /3 s is the same for all S-S edges 
in the network, regardless of their local neighborhood — the same assumption is made for the other rates 7 s , f) x , and 
7*. See also the explanation in [42] for the SIS case. 
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Figure 10: Rates for state transitions of a connected pair of node, assuming that the states of 
all their other neighbors remain unchanged. Node 1 (left-hand node) has degree k\ and has m\ 
infected enighbors that are not shown. Node 2 (right-hand node) has degree &2 an d has 777,2 
neighbors — other than node 1 — who are infected. 

equations which, along with initial conditions (|38p . can be solved numerically using standard 
methods [67] . Note that the evolution equations are completely prescribed by the functions Fk, m 
and Rk, m , and so this method can be applied to any stochastic dynamical process defined by 
transition rates Fk,m and Rk,m- For the SIS model, equations (|46l) and (|47l) were derived in [42] 
(see also [75] and [43j), with additional terms to study adaptive rewiring of the network. 

B Microscopic reversibility in equilibrium models 

We consider transition rates F^ m and Rk, m for which detailed balance holds, i.e., for which the 
dynamics exhibit microscopic reversibility [3J. Figure [TOl shows a pair of connected nodes: node 
1 — on the left — has degree k\, and node 2 — on the right — has degree k%. We assume that m\ 
neighbors of node 1, other than node 2, are in the active (infected) state; similarly, 7772 is the 
number of infected neighbors, other than node 1, of node 2. 

Starting from the state at the top of Fig. [TUl where both node 1 and node 2 are in the 
susceptible state, we consider possible cycles of state-transitions for the node pair, assuming that 
m\ and m% remain unchanged. The transition from the {S, S} state (top of figure) to the {S, 1} 
state (right of figure), for example, occurs at rate Fk 2Tn2 , since node 2 becomes infected while it 
has m2 infected neighbors. The transition from {S, 1} (right) to {/, 1} (bottom of Figure) occurs 
at rate i*]^ m i+i> since node 1 is required to become infected at a time when it has m± + 1 infected 
neighbors. The other rates are derived similarly. 

For microscopic reversibility, it is necessary that the product of the rates around a closed cycle 
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are independent of the direction of rotation around the cycle [3], I3T(, [6T] . From Fig. 1 101 this means 
we must have 

Rearranging gives the condition 

^ii,mi+l Rki.mi -^2,012 + 1 Rk2,rri2 



^iiiWi ^fci,mi+l Fk2,m,2 Rk2,rri2+1 

Since k\ and k<z (and m\ and 7712) are independent, this requires that 

Fk,m+1 Rk ,rn 



(50) 



a for all k, and all m < k, (51) 

where a is a constant. Defining u m = Fj im /Rj im , this equation can be written as the first-order 
difference equation 

u m +i = au m , (52) 

with the solution u m = a m uo. In terms of the rates, this means that detailed balance requires 
that 

Fk,m _ Fkfi 
Rk,m Rkfi 

Identifying bk as Fkp/Rkfi yields Eq. (fTTj) , 

C Derivation of equation (1201) 

Here we manipulate Eqs. (I18p and (|19j) for equilibrium models to produce the solution (|20f) and 
the algebraic equation (|2ip for p. Starting with Eq. f)19[) . we can solve for g, yielding 

' " (54) 



1 — p + pa 

Solving ([18]) for ~p k yields 



1-9 



Inserting expression (|54p into (|55p gives the desired equation (|20p . 
Next, we consider the identity (|39p . for which the PA ansatz gives 

Pk{l ~ Pk)kPk = P kPkK l ~ Ik)- (56) 

k k 

Since p k and g fc are independent of k for the case considered here — see Eq. (fl~6]) — we obtain the 
simpler expression 

p(l - 57) = 57(1 - 5) (57) 
where 57 = ( ) . Solving Eq. ([571) for 57 and using Eq. (j54"|) gives 

57 = P (1 "f + y (58) 
1 - p 2 + p 2 a 

and equating this to ^ fc ^PkPk, with Eq. (i20l) . gives Eq. (|2T1) . 
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D Critical point of spin models with up-down symmetry 



Spin models with up-down symmetry are discussed in Sec. 15.14 and it is shown there that they 
constitute a special case of the equilibrium class of models defined by relation (fT7|) . The up- 
down symmetry imposes condition (|24h on the parameters in the solution (1201) . so that /5fc can be 
expressed in the form 

(59) 



Pk 



1 



k ■ 



-p+pa / 

Recall from Sec. I5.1l that a symmetric solution of the AME with pk = 1/2 for all t exists; compar- 
ison with (|59p shows this corresponds to the solution with p = 1/(1 + \/a). Inserting this into the 
expression (|58|) for uJ — the probability that one end of a randomly chosen edge is infected — gives 
ZJ = 1/2 for the symmetric solution. 

Next, we investigate the possibility of other solutions lying near the symmetric solution; in 
the language of dynamical systems, we construct the normal form of the bifurcation [68]. We 
choose oJ as the order parameter, and let ZJ = \ + e, for small e, to explore the neighborhood of 
the symmetric solution. Equation (j58[) can then be inverted to give p in terms of e, and inserting 
this into Eq. (|2~T|) leads to a self-consistent equation for e: 



e + 



k 



1 + 



2e+ ^0 + 4(1 - a)e 2 



;i + 2e)^ 



After some rearrangement, this can be written as 



(60) 



-^-Pfctanh 



2 r z 



* l n | - + V a + 4 ( 1 



a e z 



;i + 2e)^ 



(61) 



and, in the case a = e iJ ^ T , this is very similar to Eq. (13) of [25] and Eq. (28) of [26J, which were 
derived — using very different methods (recursion method and replica trick) — for the Ising model. 

It follows that the structure of solutions of Eq. (16ip near the critical point (bifurcation point) 
can be analysed — for the class of spin models obeying (|23p — using the same methods for (|61 j) as 
used in [26, [25]. If Pk has finite fourth moment, for example, the right hand side of (I6ip can be 
expanding as a Taylor series in small e, giving 



2z 



1 



e + c 3 e 3 + C(e 5 ) as e -> 0, 



(62) 



where C3 is a coefficient involving the moments of Pk up to the fourth. Apart from the symmetric 
solution e = 0, solutions with non-zero e can be found by balancing leading-order terms, provided 
that the coefficient of e (gathering terms on the right-hand side) is positive: 



2z V y/a 



1 > 0. 



(63) 
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This inequality is equivalent to the condition 



(k 2 ) 



Thus the critical value a c is given by Eq. (|25|) . 

The case Pfe ~ A; -7 can be examined in the same fashion as in [26] 125]. The small-e expansion 
of (|6ip has leading-order terms of the form 

e ~ cie + c 3 e 3 + c e 7 " 2 as e 0. (65) 

If 7 is in the range 2 < 7 < 3, the e 7-2 term dominates both the linear and cubic terms, and 
the critical point is determined by the vanishing of coefficient c$: this occurs at a = a c = 1. For 
exponents in the range 3 < 7 < 5, the linear term dominates the e 7-2 term, and the critical point 
is again ([25|h but with a different scaling near criticality. For 7 > 5 we recover the case ([62]) . 



E PA critical point of spin models on z-regular random graphs 

For models with the symmetry ()23|) . but not necessarily possessing the equilibrium property (|17p . 
we focus here on the steady-state of the PA equations ([4]) on z-regular graphs. In particular, we 
suppose that the dynamics (through the rates i*fc !m ) depends on a parameter Q, and we derive 
the condition (128p determining the critical value (bifurcation point) of this parameter. 

We begin by noting that the property (|23p means that p = 1/2 is always a steady-state solution 
of @. Using this value, and property (|23p . on the right-hand side of the first of the PA equations 
dl]) gives the steady-state relation 

1 z 1 z 

— — ^ ^ F zz .. m Ii z ,rn (?) "I - 2 ^ ^ Fz,mBz,mip) 
m=0 m=0 

1 2 

= - g Yl l B *,™0- -9)- B z , m (p)] , (66) 

771=0 

which is satisfied — as are the other steady PA equations — if q = 1 — p in the symmetric state. 
Using this to replace q in the second of the PA equations gives, after some manipulation, the 
implicit relation for the steady-state value of p: 

j>2 f 1 - — ) Fz,mB z ,m(p) = 0. (67) 

m=0 ^ Z ' 

Next we consider the possibility of steady-states that are distinct from the symmetric state 
solution discussed above; we introduce the symbol a as a convenient shorthand for the symmetric 
state: a = {p = 1/2, p solves (I67p . q = 1 — p}. For a general — possibly non-symmetric — steady 
state, the first of the PA equations is 

z z 

= -p F Zim B z>m (l -q) + (l-p)J2 F Ztm B z>m ($). (68) 
m=0 m=0 
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We now treat p, p } and q as implicit functions of Q, the parameter denning the rates Fz,m- 
Differentiating both sides of equation (|68j) with respect to Q, evaluating at the symmetric state 
a, using the relations q\ a = 1 — p\ a and, from ([5T 



d__ 



d_ 

dQ 

dp 

dQ 



1 



P- 
—P 



dp 
dQ 



leads eventually to the relation 







} dp_ 
'dQ 



1 + 



z2p- 1 
2 1-p 



£ft 

m=0 



(69) 



(70) 



where all terms are evaluated at the symmetric state a. For this equation to be true, we must 
have either dp/dQ = or the term in square brackets must vanish (note the third term is a sum 
of positive terms so cannot be zero). The symmetry-breaking bifurcation points correspond to 
the vanishing of the square bracketed term, and this gives the critical value of p: 

z-2 



Pc 



2z 



(71) 



Finally, inserting this value into the general steady-state solution (|67p gives the criticality condi- 
tion (I28p on the rate F Z)Tn (= R ZtZ _ m ) of the PA dynamics. 



F Derivation of reduced-dimension equations for threshold mod- 
els 

Our goal here is to demonstrate that Eqs. (f30|) and (f3Tj) reduce to Eqs. (f32]) - ([M|) through an exact 
solution of (J3TT) given by (i35j) . 

We proceed to insert the ansatz (|35p into the AME (|3ip for m < (for which m values we 
have Fk )m = 0). The left hand side of (1311) then gives, after some rearrangement, 



Sk,m = (1 - Pk(0)) 



m k — m 



1 



where dots denote time derivatives. Using the identity 

1 — 6 m 



Bh , 



k,m—l 



(0) 



6 k — m + 1 
on the right hand side of Eq. (|3ip yields the condition 



-B 



k,m\ 



(72) 



(73) 



(74) 



on the function 6{t) for the ansatz (|35p to be an exact solution of (|3ip . From the initial condition 
on Sk,-m we also obtain 

0(O) = p(O) = ^-P kPk (O). (75) 
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Comparing (|32p and (|74p we see that it remains for us to show that 



f? = »m=±. (76) 

1-0 



We first prove a preliminary and rather general result, as follows. Multiplying the AME (|3ip by 
(k — m)Pk, summing over m = 0, . . . , k and over the k-classes gives: 

^X^X^ ~ m ) s Km = Pk^2 Fk,m{k - m)s Km 

km km 

-P s ^P^Yl [(* - m ) 2s k,m - (A; - m)(k - m + l)s k , m _i] . (77) 

k m 

Note that the second term on the right hand side is a telescoping series that reduces to 

-/3 S I^I^-™Km- (78) 
k m 

and that the definition of j3 s enables us to express the first term on the right hand side as 

-P S J2 Pk lL( k - m >^. (79) 
k m 

Therefore ([77]) can be solved for /3 s to yield 

-I ft Ek P k E m ( k ~ m)s k>m 



Z^k^kZ^ml^' - "V*'k,m 
1 W /„ . , \ 

(80) 



2dt 



In I ^2 P * X^ ~ m ) S k,m ) 
V k m / 



Rewriting ([Tlj) as /3 s = — A ln(l— (ft) and comparing with (|80|) gives the result that ^ k Pk E m (^~ 
m ) s k,m an d (1 — <ft) 2 are equal, up to a multiplicative constant. Using the initial conditions on 
•Sk,m and (ft determines the constant, giving 

^^k^(A:-m) Skim = z(l-</») 2 . (81) 



We now use this result and the definition of (3 s to prove (176p . For the infection probabilities 
30p . /3 s is given by 

s = Ek P kE m > A /J A; - m ) s k,m 



Ekfl<Em( fc - m )v 

Ek P k E m ( fc ~ m >k,m ~ Ek P k Em<M k ( fc ~ m ) s k,r. 
Ek P kEm( fc -"l)»k,m 

2(1 - </>) 2 " Ek ^k E m <M k (^ - m)(l - p k (0))i? fc , m 



,(1-0T • (82) 
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where the final line uses (|8ip twice and the ansatz (|35p . Dividing the numerator and denominator 
by z(l-4>) and using the identities (k-m)B k>m (cj)) = fc(l-0)B fc _i )m (0) and £ m<Afk B k-l,m{4>) = 
1 - E m >M k B k-i,m{4>) leads to 



1 



(83) 



as claimed. 

Finally, let's show that the fraction p(t) of active nodes given in the AME as p(t) = ^ k P^ Y2 m 0-' 
*k,m)) obeys equations ([32|) - ([55]) . Multiplying equation ([3"T]) by P^ and summing over m and over 
k-classes gives 



} ] Pk}, «k,m = ~ } ) Pk ^k,mSk,m - /3 S Pk ^2 ^ k ~ m ) s k,m - (fc - m + l)s kj m-l] 
km km km 

where the second term on the right hand side is easily seen to telescope to zero. Thus we have 

m = E Pk E Fk 

k m 

= E ^ k E Sk ' m ~~ E ^ k E Sk ' m 
km k m<M k 



(84) 



1-^-^^(1-^(0)) 



1- E ^(^) 



m>A/ u 



(85) 



using (|35j) . It is straightforward to verify that this reduces to equations (|32j) - (|33l . 

For completeness, note that the distribution of Sk,m for m > My. is, in general, not of the 
binomial form (|35p . which applies only to m values below the threshold M^. To obtain the values 
of 

•5k, m for these m values, note that equation (]3ip has a solution giving si^m explicitly in terms 
of 5^ m _i by using the integrating factor 



(86) 



and using (|74p for /3 s . Then we have 

Sk,m(*) = Sk,m(0) 

_ m + i) e -*(i _ 0)*-™ /" e - (i _ 0( r ))-( fc —) /3'(r) ak , TO _i(r) dr 

J o 

which can be solved recursively for m = + 1 to m = k. 



(87) 
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